



cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c
      subroutine k11(o_n,A,na,B,nb,c)
      implicit none
      integer na,nb
Cf2py intent(in) na,nb
      real*8 c
Cf2py intent(in) c
      real*8 A(na),B(nb)
Cf2py intent(in) A,B
Cf2py depend(na) A
Cf2py depend(nb) B
      real*8 o_n(na*nb)
Cf2py intent(out) o_n
Cf2py depend(na,nb) o_n
      integer i,j
      o_n=0
      do i=1,na,1
      do j=1,nb,1
      o_n(j+(i-1)*nb)=c*A(i)*B(j)
      end do
      end do
c      return
      end subroutine k11
c
c
      subroutine k21r(o_n,A,na,B,nbx,nby,c)
      implicit none
      integer na,nbx,nby
Cf2py intent(in) na,nbx,nby
      real*8 A(na),B(1:nbx,1:nby)
Cf2py intent(in) A,B
Cf2py depend(na) A
Cf2py depend(nbx,nby) B
      real*8 c
Cf2py intent(in) c
      real*8 o_n(1:nbx,1:na*nby)
Cf2py intent(out) o_n
Cf2py depend(na,nbx,nby) o_n
      integer i,i2,j2
      o_n=0
      do i=1,na,1
      do i2=1,nbx,1
      do j2=1,nby,1
      o_n(i2,j2+(i-1)*nby)=A(i)*B(i2,j2)*c
      end do
      end do
      end do
c      return
      end subroutine k21r
c
c
      subroutine k21c(o_n,A,B,c,na,nbx,nby)
      implicit none
      integer, intent(in) :: na,nbx,nby
Cf2py intent(in) na,nbx,nby
      real*8, intent(out), dimension(1:na*nbx,1:nby) :: o_n
Cf2py intent(out) o_n
Cf2py depend(na,nbx,nby) o_n
      real*8, intent(in) :: c
Cf2py intent(in) c

      real*8 A(na),B(1:nbx,1:nby)
Cf2py intent(in) A,B
Cf2py depend(na) A
Cf2py depend(nbx,nby) B
      integer i,i2,j2
      o_n=0
      do i=1,na,1
      do i2=1,nbx,1
      do j2=1,nby,1
      o_n(i2+(i-1)*nbx,j2)=A(i)*B(i2,j2)*c
      end do
      end do
      end do
c      return
      end subroutine k21c
c
c
      subroutine k22(o_n,A,nax,nay,da,B,nbx,nby,db,c)
      implicit none
      integer nax,nay,nbx,nby,da,db
Cf2py intent(in) nax,nay,nbx,nby
      real*8 A(1:nax,1:nay),B(1:nbx,1:nby)
Cf2py intent(in) A,B
Cf2py depend(nax,nay) A
Cf2py depend(nbx,nby) B
      real*8 c
Cf2py intent(in) c
      real*8 o_n(1:nax*nbx,1:nay*nby)
Cf2py intent(out) o_n
Cf2py depend(nax,nay,nbx,nby) o_n
      integer i1,i2,j1,j2,io,jo
      o_n=0

      do i1=1,nax,1
        if (da.eq.0) then
          do j1=1,nay,1
      do i2=1,nbx,1
        if (db.eq.0) then
          do j2=1,nby,1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end do
        else if (db.eq.1) then
          j2=i2
          io=i2+(i1-1)*nbx
          jo=j2+(j1-1)*nby
          o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
        else if (db.eq.2) then
          if (i2.eq.1) then
            j2=2
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else if (i2.eq.nbx) then
            j2=nbx-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else
            j2=i2+1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
            j2=i2-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end if
        end if
      end do
          end do
        else if (da.eq.1) then
          j1=i1
      do i2=1,nbx,1
        if (db.eq.0) then
          do j2=1,nby,1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end do
        else if (db.eq.1) then
          j2=i2
          io=i2+(i1-1)*nbx
          jo=j2+(j1-1)*nby
          o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
        else if (db.eq.2) then
          if (i2.eq.1) then
            j2=2
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else if (i2.eq.nbx) then
            j2=nbx-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else
            j2=i2+1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
            j2=i2-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end if
        end if
      end do
        else if (da.eq.2) then
          if (i1.eq.1) then
            j1=2
      do i2=1,nbx,1
        if (db.eq.0) then
          do j2=1,nby,1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end do
        else if (db.eq.1) then
          j2=i2
          io=i2+(i1-1)*nbx
          jo=j2+(j1-1)*nby
          o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
        else if (db.eq.2) then
          if (i2.eq.1) then
            j2=2
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else if (i2.eq.nbx) then
            j2=nbx-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else
            j2=i2+1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
            j2=i2-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end if
        end if
      end do
          else if (i1.eq.nax) then
            j1=nax-1
      do i2=1,nbx,1
        if (db.eq.0) then
          do j2=1,nby,1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end do
        else if (db.eq.1) then
          j2=i2
          io=i2+(i1-1)*nbx
          jo=j2+(j1-1)*nby
          o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
        else if (db.eq.2) then
          if (i2.eq.1) then
            j2=2
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else if (i2.eq.nbx) then
            j2=nbx-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else
            j2=i2+1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
            j2=i2-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end if
        end if
      end do
          else
            j1=i1+1
      do i2=1,nbx,1
        if (db.eq.0) then
          do j2=1,nby,1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end do
        else if (db.eq.1) then
          j2=i2
          io=i2+(i1-1)*nbx
          jo=j2+(j1-1)*nby
          o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
        else if (db.eq.2) then
          if (i2.eq.1) then
            j2=2
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else if (i2.eq.nbx) then
            j2=nbx-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else
            j2=i2+1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
            j2=i2-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end if
        end if
      end do
            j1=i1-1
      do i2=1,nbx,1
        if (db.eq.0) then
          do j2=1,nby,1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end do
        else if (db.eq.1) then
          j2=i2
          io=i2+(i1-1)*nbx
          jo=j2+(j1-1)*nby
          o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
        else if (db.eq.2) then
          if (i2.eq.1) then
            j2=2
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else if (i2.eq.nbx) then
            j2=nbx-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          else
            j2=i2+1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
            j2=i2-1
            io=i2+(i1-1)*nbx
            jo=j2+(j1-1)*nby
            o_n(io,jo)=A(i1,j1)*B(i2,j2)*c
          end if
        end if
      end do
          end if
        end if
      end do

c      return
      end subroutine k22
c
c
c
c
c

      subroutine kk22(o_n,na,nb,nc,da,db,dc,A,B,C,g)
c kron(a,kron(b,c)
      implicit none
      integer na,nb,nc,da,db,dc
Cf2py intent(in) na,nb,nc,da,db,dc
      real*8 A(1:na,1:na),B(1:nb,1:nb),C(1:nc,1:nc)
Cf2py intent(in) A,B,C
Cf2py depend(na) A
Cf2py depend(nb) B
Cf2py depend(nc) C
      real*8 g
Cf2py intent(in) g
      real*8 o_n(1:na*nb*nc,1:na*nb*nc)
Cf2py intent(out) o_n
Cf2py depend(na,nb,nc) o_n
      integer ia,ib,ic,ja,jb,jc,io,jo
      o_n=0
c      do ia=1,na,1
c      do ib=1,nb,1
c      do ic=1,nc,1
c      do ja=1,na,1
c      do jb=1,nb,1
c      do jc=1,nc,1
c      io=ic+(ib-1)*nc+(ia-1)*nb*nc
c      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
c      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
c      end do
c      end do
c      end do
c      end do
c      end do
c      end do

      do ia=1,na,1
        if (da.eq.0) then
          do ja=1,na,1
      do ib=1,nb,1
        if (db.eq.0) then
          do jb=1,nb,1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end do
        else if (db.eq.1) then
          jb=ib

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
        else if (db.eq.2) then
          if (ib.eq.1) then
            jb=2

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else if (ib.eq.nb) then
            jb=nb-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else
            jb=ib+1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
            jb=ib-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end if
        end if
      end do
          end do
        else if (da.eq.1) then
          ja=ia
      do ib=1,nb,1
        if (db.eq.0) then
          do jb=1,nb,1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end do
        else if (db.eq.1) then
          jb=ib

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
        else if (db.eq.2) then
          if (ib.eq.1) then
            jb=2

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else if (ib.eq.nb) then
            jb=nb-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else
            jb=ib+1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
            jb=ib-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end if
        end if
      end do
        else if (da.eq.2) then
          if (ia.eq.1) then
            ja=2
      do ib=1,nb,1
        if (db.eq.0) then
          do jb=1,nb,1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end do
        else if (db.eq.1) then
          jb=ib

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
        else if (db.eq.2) then
          if (ib.eq.1) then
            jb=2

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else if (ib.eq.nb) then
            jb=nb-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else
            jb=ib+1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
            jb=ib-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end if
        end if
      end do
          else if (ia.eq.na) then
            ja=na-1
      do ib=1,nb,1
        if (db.eq.0) then
          do jb=1,nb,1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end do
        else if (db.eq.1) then
          jb=ib

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
        else if (db.eq.2) then
          if (ib.eq.1) then
            jb=2

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else if (ib.eq.nb) then
            jb=nb-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else
            jb=ib+1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
            jb=ib-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end if
        end if
      end do
          else
            ja=ia+1
      do ib=1,nb,1
        if (db.eq.0) then
          do jb=1,nb,1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end do
        else if (db.eq.1) then
          jb=ib

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
        else if (db.eq.2) then
          if (ib.eq.1) then
            jb=2

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else if (ib.eq.nb) then
            jb=nb-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else
            jb=ib+1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
            jb=ib-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end if
        end if
      end do
            ja=ia-1
      do ib=1,nb,1
        if (db.eq.0) then
          do jb=1,nb,1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end do
        else if (db.eq.1) then
          jb=ib

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
        else if (db.eq.2) then
          if (ib.eq.1) then
            jb=2

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else if (ib.eq.nb) then
            jb=nb-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          else
            jb=ib+1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
            jb=ib-1

      do ic=1,nc,1
        if (dc.eq.0) then
          do jc=1,nc,1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end do
        else if (dc.eq.1) then
          jc=ic
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
        else if (dc.eq.2) then
          if (ic.eq.1) then
            jc=2
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else if (ic.eq.nc) then
            jc=nc-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          else
            jc=ic+1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
            jc=ic-1
      io=ic+(ib-1)*nc+(ia-1)*nb*nc
      jo=jc+(jb-1)*nc+(ja-1)*nb*nc
      o_n(io,jo)=A(ia,ja)*B(ib,jb)*C(ic,jc)*g
          end if
        end if
      end do
          end if
        end if
      end do
          end if
        end if
      end do
c      return
      end subroutine kk22











c
c
      subroutine sk22(o_n,na,nb,ca,cb,da,db,dia,dib,A,B,IA,IB,g)
c o=kron22(A,IB)
c o+=kron22(IA,B)
      implicit none
      integer na,nb,da,db,dia,dib
Cf2py intent(in) na,nb
      real*8 A(1:na,1:na),B(1:nb,1:nb)
      real*8 IA(1:na,1:na),IB(1:nb,1:nb)
Cf2py intent(in) A,B,IA,IB
Cf2py depend(na) A,IA
Cf2py depend(nb) B,IB
      real*8 ca,cb,g
Cf2py intent(in) ca,cb,g
      real*8 o_inter(1:na*nb,1:na*nb), o_n(1:na*nb,1:na*nb)
Cf2py intent(out) o_n
Cf2py depend(na,nb) o_n,o_inter
      o_n=0
      o_inter=0
      call k22(o_inter,A,na,na,da,IB,nb,nb,dib,ca)
      call k22(o_n,IA,na,na,dia,B,nb,nb,db,cb)
      o_n=(o_n+o_inter)*g
c      return
      end subroutine sk22
c
c
      subroutine sskk22(o_n,na,nb,nc,ca,cb,cc,d,A,B,C,IA,IB,IC,g)
c o=kk22(A,kk22(IB,IC))
c o+=kk22(IA,kk22(B,IC))
c +=kk22(IA,kk22(IB,C))
c kron(a,kron(b,c)
      implicit none
      integer na,nb,nc,d(6)
Cf2py intent(in) na,nb,nc,d
      real*8 A(1:na,1:na),B(1:nb,1:nb),C(1:nc,1:nc)
      real*8 IA(1:na,1:na),IB(1:nb,1:nb),IC(1:nc,1:nc)
Cf2py intent(in) A,B,C,IA,IB,IC
Cf2py depend(na) A,IA
Cf2py depend(nb) B,IB
Cf2py depend(nc) C,IC
      real*8 ca,cb,cc,g
Cf2py intent(in) ca,cb,cc,g
      real*8 o_n(1:na*nb*nc,1:na*nb*nc)
      real*8 o_inter1(1:na*nb*nc,1:na*nb*nc)
      real*8 o_inter2(1:na*nb*nc,1:na*nb*nc)
Cf2py intent(out) o_n
Cf2py depend(na,nb,nc) o_n
      o_n=0
      o_inter1=0
      o_inter2=0
      call kk22(o_inter1,na,nb,nc,d(1),d(5),d(6),A,IB,IC,ca)
      call kk22(o_inter2,na,nb,nc,d(4),d(2),d(6),IA,B,IC,cb)
      call kk22(o_n,na,nb,nc,d(4),d(5),d(3),IA,IB,C,cc)
      o_n=g*(o_n+o_inter1+o_inter2)
c      return
      end subroutine sskk22
c
c
      subroutine sp(o_n,na,nb,cx,cy,cz,d,Ax,Ay,Az,Bx,By,Bz,c)
      implicit none
c scalar product vector of n component of tensor rank da db
c declaration of dimensions
      integer na,nb,d(6)
Cf2py intent(in) na,nb,d
      real*8 c,cx,cy,cz
Cf2py intent(in) c,cx,cy,cz
      real*8 o_n(1:na*nb,1:na*nb)
      real*8 o_inter1(1:na*nb,1:na*nb),o_inter2(1:na*nb,1:na*nb)
Cf2py intent(out) o_n
Cf2py depend(na,nb) o_n
      real*8 Ax(1:na,1:na),Ay(1:na,1:na),Az(1:na,1:na)
Cf2py intent(in) Ax,Ay,Az
Cf2py depend(na) Ax,Ay,Az
      real*8 Bx(1:nb,1:nb),By(1:nb,1:nb),Bz(1:nb,1:nb)
Cf2py intent(in) Bx,By,Bz
Cf2py depend(nb) Bx,By,Bz
      o_inter1=0
      o_inter2=0
      o_n=0
      call k22(o_inter1,Ax,na,na,d(1),Bx,nb,nb,d(4),cx)
      call k22(o_inter2,Ay,na,na,d(2),By,nb,nb,d(5),cy)
      call k22(o_n,Az,na,na,d(3),Bz,nb,nb,d(6),cz)
      o_n=c*(o_n+o_inter1+o_inter2)
c      return
      end subroutine sp

      subroutine sptv(o_n,na,Ax,Ay,Az,B,c)
      implicit none
c scalar product vector of n component of tensor rank da db
c declaration of dimensions
      integer na
Cf2py intent(in) da
      real*8 o_n(1:na,1:na)
Cf2py intent(out) o_n
Cf2py depend(na) o_n
      real*8 Ax(1:na,1:na),Ay(1:na,1:na),Az(1:na,1:na)
Cf2py intent(in) Ax,Ay,Az
Cf2py depend(na) Ax,Ay,Az
      real*8 B(3)
      real*8 c
Cf2py intent(in) B,c
      o_n=0
      o_n=B(1)*Ax
      o_n=o_n+B(2)*Ay
      o_n=o_n+B(3)*Az
      o_n=o_n*c
c      return
      end subroutine sptv
